Qual é a alturas média das årvores em uma floresta?
Departamento de Ecologia IB-USP
Qual é a alturas média das årvores em uma floresta?
Fazer afirmaçÔes sobre as caracterĂsticas de uma população, com base em informaçÔes dadas por amostras.
População Ă© qualquer conjunto de elementos que vocĂȘ estĂĄ observando
Amostra é qualquer subconjunto desta população.
## [1] 22.72789
## [1] 32.3065
## [1] 5.683881
A variĂĄvel aleatĂłria X Ă© uma variĂĄvel que tem um valor Ășnico (determinado aleatoriamente) para cada resultado de um experimento (lembre-se esse termo veio da teoria de probabilidade).
Exemplos:
Altura das ĂĄrvores de uma floresta (!!)
NĂșmero de presas capturadas por um predador em um determinado dia;
Comprimento de um peixe adulto selecionado aleatoriamente.
As variĂĄveis aleatĂłrias podem ser discretas ou contĂnuas.
NĂșmero ou a quantidade observada na unidade experimental ou tentativa.
Representada por nĂșmeros inteiros (0, 1, 2, 3, 4…);
NĂŁo pode conter nĂșmeros negativos;
NĂșmero finito de possibilidades;
Podemos achar a probabilidade de cada evento.
Ex: peso, altura, distĂąncia, pH, biomassa, etc.
Representada por nĂșmeros nĂŁo inteiros (1.3, - 1.54, - 1.7);
Pode conter nĂșmeros negativos;
NĂșmero infinito de possibilidades;
Probabilidade de cada evento Ă© zero.
Associa cada possĂvel valor da variĂĄvel (X) Ă sua probabilidade de ocorrĂȘncia P(X).
Quando conhecemos todos os valores de uma variåvel aleatória, juntamente com suas respectivas probabilidades, temos uma distribuição de probabilidades.
Jogar moedas: cara 0, coroa 1
Presença/ ausĂȘncia de uma espĂ©cie em locais amostrados: 0 - 1
ParĂąmetro: p - probabilidade de sucesso
E[X] = p
var[X] = p(1-p)
E[X] = np
e a variĂąncia Ă©
var[X] = np(1-p)
NĂșmero sementes predadas em cada experimento contendo 20 semente.
NĂșmero de animais infectados em relação ao nĂșmero de capturados em cada local
Probabilidade de predação de girinos por odonatas. P = 0.30
Temos 6 girinos num lago.
Qual robabilidades de que 0, 1, 2, 3, 5 ou todos sejam predados?
Distribuição de probabilidade discreta.
Probabilidade de uma sĂ©rie de eventos ocorrem em um perĂodo fixo de tempo, ĂĄrea, volume, quadrante, etc.
E[X] = (λ)
var[X] = (λ)
Abundùncia de espécies em um local
NĂșmero de indivĂduos capturados por tempo
â Metade (50%) estĂĄ acima (e abaixo) da mĂ©dia
â Aproximadamente 68% estĂĄ dentro de 1 desvio padrĂŁo da mĂ©dia
â Aproximadamente 95% estĂĄ dentro de 2 desvios padrĂ”es da mĂ©dia
â Virtualmente todos os valores estĂŁo dentro de 3 desvios padrĂ”es da mĂ©dia
"toda soma de variĂĄveis aleatĂłrias independentes de mĂ©dia finita e variĂąncia limitada Ă© aproximadamente Normal, desde que o nĂșmero de termos da soma seja suficientemente grande".
Independentemente do tipo de distribuição da população, na medida em que o tamanho da amostra aumenta, a distribuição das médias amostrais tende a uma distribuição Normal.
Qual Ă© a probabilidade de que um peixe capturado aleatoriamente tenha 20,15 cm ou mais?
Média da população é 17,1 cm e o desvio padrão é de 1,21 cm?
prob <- pnorm(q = 20.15, mean = 17.1, sd = 1.21, lower.tail = F) prob
## [1] 0.005856729
levels(arv$sp)
## [1] "sp1" "sp2"
Serå que a altura média de uma espécie nessa floresta é diferente da altura média da outra espécie?
Ou seja, serå que as espécies pertencem a uma mesma distribuição e as diferenças encontradas são devido ao acaso ou serå que cada espécie é uma variåvel aleatória com médias diferentes?
library(lattice)
histogram(~alt|sp,
data = arv,
layout = c(1,2),
strip=F,
strip.left=T)
mean(arv$alt[arv$sp=="sp2"]) - mean(arv$alt[arv$sp=="sp1"])
## [1] 5.109466
Posso afirmar que as médias são diferentes?
H0 - não hå diferença entre as alturas das duas espécies. A média da altura da espécie 1 é igual à média da altura da espécie 2 e a diferença observada se deve ao acaso.
H1 - hå diferença entre as alturas das duas espécies. A média da altura da espécie 1 é diferente da média da altura da espécie 2.
AFIRMAĂĂO: nĂŁo hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: hå diferenças.
AFIRMAĂĂO: hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: não hå diferenças.
ERRO DO TIPO I (α): rejeitar a hipótese nula dada que ela é verdadeira.
ERRO DO TIPO II: aceitar a hipĂłtese nula dado que ela Ă© falsa.
Contrastar a probabilidade de que a diferença que vocĂȘ encontrou nas mĂ©dias das alturas das duas espĂ©cies pode se considerada devido ao acaso ou nĂŁo.
Qual é a probabilidade de que as duas amostras pertencem à mesma população de medidas?
Podemos ou não rejeitar a hipótese nula (de que a diferença é ao acaso), baseado em alguma probabilidade de estarmos cometendo algum erro (α).
"O nĂșmero de observaçÔes independentes menos o nĂșmero de parĂąmetros estimados."
Quantas observaçÔes independentes foram usadas para calcular a média e a variùncia dos dados sob a hipótese nula?
100 ĂĄrvores
2 parĂąmetros
98 graus de liberdade
NĂŁo entraremos em detalhes de como se calcula a estatĂstica t e como Ă© a sua distribuição de probabilidades (recomendo ver na wikipedia). O importante a saber aqui Ă© que tendo as duas amostras modeladas como uma distribuição normal, calculamos a estatĂstica t e com um certo valor de graus de liberdade, nĂłs podemos recorrer Ă s tabelas da distribuição t para ver qual a probabilidade de se obter aquele valor que obtivemos dada que a hipĂłtese nula Ă© verdadeira. Se esta probabilidade for muito pequena, a chance de estarmos caindo no erro do tipo 1 Ă© pequena.
No R, a função usada para fazer um teste t é t.test.
t.test(arv$alt~arv$sp, var.equal=T)
## ## Two Sample t-test ## ## data: arv$alt by arv$sp ## t = -5.0125, df = 98, p-value = 2.387e-06 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -7.132311 -3.086621 ## sample estimates: ## mean in group sp1 mean in group sp2 ## 20.17316 25.28262
#o argumento var.equa=T diz que as variùncias das duas espécies são as mesmas # podemos fazer o teste com a premissa de que as variùncias não sejam iguais, # usando var.equal=F, isso vai "consumir" graus de liberdade.
Vamos interpretar este resultado: temos o valor de t calculado e os graus de liberade (jå calculado pela gente antes). Com estes 2 valores, o que o R faz é calcular a probabilidade de que a hipótese nula seja verdadeira. Esse é o valor-p, que neste exemplo deu um valor beem pequeno. Vemos então que podemos rejeitar H0, com uma probabilidade muito baixa de estarmos incorrendo no erro de tipo 1, e desta forma acreditamos que as espécies de årvores tenham distribuiçÔes de alturas diferentes (médias diferentes).
E se tivĂ©ssemos na verdade 3 espĂ©cies de ĂĄrvores para compararmos as alturas? Acabamos de ver que o teste t serve apenas para comparar duas amostras. O que chamamos de modelos lineares Ă© uma grande famĂlia de modelos estatĂsticos que modelam uma relação linear entre a variĂĄvel dependente (Y) e a(s) variĂĄvel(is) independente(s). Nesta famĂlia encontra-se o teste t que acabamos de ver, a AnĂĄlise de VariĂąnci (ANOVA), a regressĂŁo, etc.
Chamamos de Modelos Lineare Gerais (LMs), aqueles que possuem como premissa a variĂĄvel dependente Y como de distribuição Normal, e de Modelos Lineares Genearlizados (GLMs) os que assumem outros tipos de distribuição para a variĂĄvel dependente (como a Binomial e Poisson). No R, a função que usamos para estes modelos sĂŁo lm e glm (hĂĄ alguns mais especĂficos como o prĂłrpio t.test e o aov para ANOVA). O R lida com todos os modelos lineares de maneira similar, e um detalhe importante Ă© que todas as anĂĄlises desta famĂlia devem ser salvas em um objeto, para seu resultado ser apresentado apĂłs o comando summary.
Ă importante lembrar que neste roteiro estamos lidando com modelos paramĂ©tricos, ou seja, que modela a variĂĄvel de interesse como uma variĂĄvel aleatĂłria pertencente a uma certa distribuição de probabilidades. Existem outras abordagens estatĂsticas nĂŁo-paramĂ©tricas que fazem testes que nĂŁo tem como premissa a distribuição de probabilidades. Veja um exemplo de teste t nĂŁo paramĂ©trico neste roteiro.
Mesmo fazendo parte de uma mesma famĂlia de modelos, nas sessĂ”es seguintes vamos separar os modelos lineares pelos nomes tradicionais das anĂĄlises: anĂĄlise de variĂąncia, regressĂŁo linear simples e anĂĄlise de covariĂąncia. Depois, veremos alguns exemplos de GLMs aplicĂĄveis a dados ecolĂłgicos.
Vamos entĂŁo incluir no nosso exemplo inicial uma terceira espĂ©cie amostrada (e vamos diminuir o nĂșmero de amostras por espĂ©cie para ficar mais fĂĄcil de fazer os cĂĄlculos passo a passo):
alt.sp3 <- c(36.90076761,33.23131727,32.82767311,24.93410577,31.87626329,
28.76824248,26.6436144,31.41238398,27.65383929,31.70425719,
33.60752445,28.30263354,29.32210674,23.4790054,22.20793046,
18.28797293,23.39044736,28.75820917,29.94344703,34.47430235,
32.09820862,26.99845592,27.64858951,31.90557306,33.63236368,
25.35966179,35.00885172,32.89125598,18.98802229,30.8516906,
22.04191398,35.37232889,35.67150408,34.30257037,24.43038475,
20.6274663,31.15717916,23.27023231,36.20568545,34.33068448,
28.40769054,31.00578874,29.50989158,31.69915991,37.25975797,
30.1279997,38.623371,35.36472829,29.88304259,31.00145491)
sp3 <- data.frame(alt=alt.sp3, sp=rep("sp3",50))
arv2 <- rbind(arv,sp3)
amostras <- seq(1,150,by=5)
arv3 <- arv2[ amostras, ]
Neste exemplo lidaremos com um modelo simples de anova chamado one way ANOVA, porque lida apenas com uma variĂĄvel independente categĂłrica (fator).
A ANOVA testarå a hipótese nula de que as médias das alturas das 3 espécies não diferem.
Faremos agora a ANOVA "na unha", para entendermos melhor cada passo da anålise e a interpretação de seus resultados para depois usarmos a função do R que faz esta anålise diretamente.
Vamos calcular os valores da tabela de ANOVA. Começando com os desvios quadrĂĄticos, ou seja, quanto os dados desviam da mĂ©dia (idĂ©ia parecida com a variĂąncia). O ponto importante Ă© que essa variação Ă© aditiva e portanto, pode se decomposta. A variação total Ă© decomposta em variação relacionada ao tratamento (entre grupos), no nosso caso Ă s espĂ©cies, e uma variação interna dos grupos (chamada de erro). A estatĂstica F Ă© a razĂŁo entre essas variaĂ”es apĂłs dividir cada uma delas pelos seus respectivos graus de liberdade. Complicou? Vamos fazer os cĂĄlculos e ver os grĂĄficos para ver se entendemos melhor:
A tabela de ANOVA que vamos construir Ă© essa:
Primeiro vamos mudar a forma do nosso data frame para facilitar os comandos
arv4 <- data.frame(sp1=arv3$alt[arv3$sp=="sp1"],
sp2=arv3$alt[arv3$sp=="sp2"],
sp3=arv3$alt[arv3$sp=="sp3"])
arv4
## sp1 sp2 sp3 ## 1 21.282000 19.14508 36.90077 ## 2 19.875790 18.66252 28.76824 ## 3 23.562570 25.10581 33.60752 ## 4 19.313744 18.10266 18.28797 ## 5 26.308616 21.48084 32.09821 ## 6 8.301795 28.08559 25.35966 ## 7 27.192502 20.47805 22.04191 ## 8 15.941866 29.22927 20.62747 ## 9 20.596624 19.32911 28.40769 ## 10 22.048945 18.21446 30.12800
boxplot(arv4)
Vamos calcular a média geral, as médias para cada espécie, e as diferenças entre a média geral e a média para cada espécie, para chegar ao valor da soma dos quadrados total:
media.geral <- mean(arv3$alt) media.geral
## [1] 23.28284
medias.sps <- apply(arv4, 2, mean) medias.sps
## sp1 sp2 sp3 ## 20.44245 21.78334 27.62274
dif.geral <- arv4 - media.geral ss.especies <- dif.geral^2 ss.especies
## sp1 sp2 sp3 ## [1,] 4.00337202 17.121094 185.447876 ## [2,] 11.60801175 21.347420 30.089610 ## [3,] 0.07824755 3.323227 106.599051 ## [4,] 15.75374565 26.834280 24.948725 ## [5,] 9.15530094 3.247214 77.710675 ## [6,] 224.43179074 23.066347 4.313177 ## [7,] 15.28543337 7.866884 1.539904 ## [8,] 53.88994194 35.359990 7.051024 ## [9,] 7.21577099 15.631997 26.264064 ## [10,] 1.52250306 25.688498 46.856173
ss.total <- sum(ss.especies) ss.total
## [1] 1033.251
Para calcular os desvios quadrĂĄticos totais, nĂłs subtraĂmos cada altura das ĂĄrvores pela mĂ©dia geral e elevamos ao quadrado. Veja o grĂĄfico:
vetor.cor <- rep(1:3, each=10) #vetor de cores
plot(x = c(1:30), y = arv3$alt, ylim=c(10,40),
pch=(rep(c(15,16,17), each=10)),
col=vetor.cor,
ylab="Variåvel Resposta", xlab="ObservaçÔes",
main="Variação total")
for(i in 1:30)
{
lines(c(i,i),c(arv3$alt[i],mean(arv3$alt)),col=vetor.cor[i])
}
abline(h=media.geral)
Agora vamos fazer a somatĂłria dos desvios quadrĂĄticos dentro de cada grupo (ss.intra).
O grĂĄfico para entender esse cĂĄlculo:
vetor.medias<-rep(medias.sps, each=10)
plot(c(1:30), arv3$alt, ylim=c(10,40),
pch=(rep(c(15,16,17),each=10)),
col=vetor.cor,
main="Variação Intra Grupos",
ylab="Variåvel Resposta", xlab="ObservaçÔes")
for(i in 1:30)
{
lines(c(i,i),c(vetor.medias[i],arv3$alt[i]),col=vetor.cor[i])
}
lines(c(1,10),c(medias.sps[1],medias.sps[1]),col=1)
lines(c(11,20),c(medias.sps[2],medias.sps[2]),col=2)
lines(c(21,30),c(medias.sps[3],medias.sps[3]),col=3)
CĂĄlculo dos valores:
medias.sps
## sp1 sp2 sp3 ## 20.44245 21.78334 27.62274
ss.sp1=sum((arv4$sp1-medias.sps["sp1"])^2) ss.sp1
## [1] 262.2655
ss.sp2=sum((arv4$sp2-medias.sps["sp2"])^2) ss.sp2
## [1] 157.0018
ss.sp3=sum((arv4$sp3-medias.sps["sp3"])^2) ss.sp3
## [1] 322.4728
ss.intra=ss.sp1+ss.sp2+ss.sp3 ss.intra
## [1] 741.7401
A soma dos quadrados entre grupos:
plot(c(1:30), vetor.medias, ylim=c(10,40),
pch=(rep(c(15,16,17),each=10)),
col=vetor.cor,
main="Variação Entre Grupos",
ylab="Variåvel Resposta", xlab="ObservaçÔes")
for(i in 1:30)
{
lines(c(i,i),c(vetor.medias[i],mean(vetor.medias)),col=vetor.cor[i])
}
abline(h=media.geral)
points(c(1:30),arv3$alt, ylim=c(10,50),
pch=(rep(c(0,1,2),each=10)), col=vetor.cor, cex=0.5)
medias.sps
## sp1 sp2 sp3 ## 20.44245 21.78334 27.62274
media.geral
## [1] 23.28284
ss.entre=10*sum((medias.sps-media.geral)^2) ss.entre
## [1] 291.5112
#conferindo os cĂĄlculos ss.intra+ss.entre
## [1] 1033.251
ss.total
## [1] 1033.251
CĂĄlculo do F
# Desvios médios ms.entre=ss.entre/2 ms.intra=ss.intra/27 ms.entre
## [1] 145.7556
ms.intra
## [1] 27.47186
#F - razĂŁo das variĂąncias F.sps=ms.entre/ms.intra F.sps
## [1] 5.305634
Vamos ver na distribuição F qual a probabilidade de termos econtrado o valor F.sps ao acaso:
curve(expr=df(x, 2,27),
main="Distribuição F de Fisher (df=2,27)",
xlab="Valor F",
ylab="Densidade ProbabilĂstica (df)",
xlim=c(2,12))
abline(v=F.sps, col="red")
abline(h=0, lty=2)
xf=seq(F.sps, 12, 0.01)
ydf=df(xf, 2, 27)
polygon(c(F.sps,xf),c(0,ydf),col="red")
text(x= 7,y=0.08,paste("pf(x) =",
round(pf(F.sps,2,27,lower.tail=F),4)),
cex=0.8, col="red")
CĂĄlculo do P
p.sps=pf(F.sps, 2, 27, lower.tail=FALSE) p.sps
## [1] 0.01139248
Para mais informaçÔes sobre o F e as comparaçÔes com o teste t, veja esse roteiro.
Colocando os dados calculados em nossa tabela de anova
A função que constrĂłi esta tabela de ANOVA e faz o teste estatĂstico Ă© a aov. Vamos comparar:
anova.sps <- aov(alt~sp, data=arv3) summary(anova.sps)
## Df Sum Sq Mean Sq F value Pr(>F) ## sp 2 291.5 145.76 5.306 0.0114 * ## Residuals 27 741.7 27.47 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Agora nĂłs sabemos de ondem vem os nĂșmeros nesta tabela e como foi calculado o P, certo?
Onde estão as diferenças?
veja os grĂĄficos
faça teste à posteriori - Tukey HDS
TukeyHSD(anova.sps)
## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = alt ~ sp, data = arv3) ## ## $sp ## diff lwr upr p adj ## sp2-sp1 1.340893 -4.4708805 7.152667 0.8360276 ## sp3-sp1 7.180300 1.3685259 12.992073 0.0132022 ## sp3-sp2 5.839406 0.0276327 11.651180 0.0487468
#revendo grĂĄfico dos dados
stripchart( arv3$alt~arv3$sp,
vertical = TRUE, pch = 16,
col = c("black", "red", "green"))
#boxplot boxplot(arv4)
lm.anova.sp <- lm(alt~sp, data=arv3) summary.aov(lm.anova.sp) #mesmo que o summary do aov
## Df Sum Sq Mean Sq F value Pr(>F) ## sp 2 291.5 145.76 5.306 0.0114 * ## Residuals 27 741.7 27.47 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mostra o resultado de forma um pouco diferente, consegue entender? summary(lm.anova.sp)
## ## Call: ## lm(formula = alt ~ sp, data = arv3) ## ## Residuals: ## Min 1Q Median 3Q Max ## -12.1407 -3.0002 -0.0742 3.2719 9.2780 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 20.442 1.657 12.334 1.32e-12 *** ## spsp2 1.341 2.344 0.572 0.57202 ## spsp3 7.180 2.344 3.063 0.00492 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 5.241 on 27 degrees of freedom ## Multiple R-squared: 0.2821, Adjusted R-squared: 0.229 ## F-statistic: 5.306 on 2 and 27 DF, p-value: 0.01139
plot(arv3$alt ~ nutri)
alt.nutr <- lm(alt ~ nutri, data = arv3) summary(alt.nutr)
## ## Call: ## lm(formula = alt ~ nutri, data = arv3) ## ## Residuals: ## Min 1Q Median 3Q Max ## -6.7176 -1.7171 0.0088 1.4124 9.2072 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 14.38711 1.05832 13.594 7.44e-14 *** ## nutri 0.58006 0.05971 9.714 1.82e-10 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.906 on 28 degrees of freedom ## Multiple R-squared: 0.7712, Adjusted R-squared: 0.763 ## F-statistic: 94.37 on 1 and 28 DF, p-value: 1.818e-10
plot(arv3$alt ~ nutri) abline(alt.nutr, col = "red") # reta da regressĂŁo #reta da hipĂłtese nula de ausĂȘncia de efeito dos nutrientes abline(h = mean(arv3$alt), col = "blue", lty = 2)
Y continuo
X1 contĂnuo
X2 categĂłrico
plot(alt~nutri, data=arv3, col=arv3$sp)
alt.sp.nutr <- lm(alt ~ nutri + sp, data = arv3) summary(alt.sp.nutr)
## ## Call: ## lm(formula = alt ~ nutri + sp, data = arv3) ## ## Residuals: ## Min 1Q Median 3Q Max ## -5.6682 -0.9554 -0.1585 1.3581 7.6134 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 13.3971 1.1034 12.142 3.23e-12 *** ## nutri 0.5256 0.0561 9.369 8.09e-10 *** ## spsp2 1.6421 1.1423 1.437 0.16251 ## spsp3 3.8326 1.1965 3.203 0.00357 ** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.553 on 26 degrees of freedom ## Multiple R-squared: 0.836, Adjusted R-squared: 0.817 ## F-statistic: 44.16 on 3 and 26 DF, p-value: 2.401e-10
alt.sp.nutr2 <- lm(alt ~ nutri + sp + nutri:sp, data=arv3) alt.sp.nutr3 <- lm(alt ~ nutri*sp, data=arv3) summary(alt.sp.nutr3)
## ## Call: ## lm(formula = alt ~ nutri * sp, data = arv3) ## ## Residuals: ## Min 1Q Median 3Q Max ## -5.6724 -0.9444 -0.3313 1.0416 7.0702 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 13.40167 1.42727 9.390 1.66e-09 *** ## nutri 0.52527 0.08906 5.898 4.38e-06 *** ## spsp2 2.93717 1.97080 1.490 0.149 ## spsp3 0.43636 2.75666 0.158 0.876 ## nutri:spsp2 -0.10095 0.12423 -0.813 0.424 ## nutri:spsp3 0.17187 0.14350 1.198 0.243 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 2.474 on 24 degrees of freedom ## Multiple R-squared: 0.8578, Adjusted R-squared: 0.8282 ## F-statistic: 28.96 on 5 and 24 DF, p-value: 2.003e-09
summary.aov(alt.sp.nutr3)
## Df Sum Sq Mean Sq F value Pr(>F) ## nutri 1 796.8 796.8 130.179 3.53e-11 *** ## sp 2 66.9 33.5 5.467 0.0111 * ## nutri:sp 2 22.6 11.3 1.846 0.1796 ## Residuals 24 146.9 6.1 ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Outras distribuiçÔes de probabilidades além da Normal
VariĂąncias nĂŁo homogĂȘneas
Mais comuns:
Poisson â Ășteis para dados de contagem
Binomial â Ășteis para dados com proporçÔes, ou presença/ausĂȘncia
Uma hipótese sobre a distribuição da variåvel resposta Yi.
Especificação do modelo (variåveis X).
3.Função de ligação entre a média Yi e as variåveis X. Função que lineariza a regressão
| Distribuição | link function |
|---|---|
| normal | identity |
| poisson | log |
| binomial | logit |
| Gamma | reciprocal |
Y - NĂșmero de atropelamentos de anuros em uma estrada
X - distĂąncia a um Parque Nacional.
plot( atrop$TOT.N ~ atrop$D.PARK,
xlab = " DistĂąncia ao parque",
ylab = "NĂșmero de atropelamentos")
mod <- glm(TOT.N ~ D.PARK, data = atrop,
family = poisson)
summary(mod)
## ## Call: ## glm(formula = TOT.N ~ D.PARK, family = poisson, data = atrop) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -8.1100 -1.6950 -0.4708 1.4206 7.3337 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 4.316e+00 4.322e-02 99.87 <2e-16 *** ## D.PARK -1.059e-04 4.387e-06 -24.13 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 1071.4 on 51 degrees of freedom ## Residual deviance: 390.9 on 50 degrees of freedom ## AIC: 634.29 ## ## Number of Fisher Scoring iterations: 4
anova(mod, test = "Chisq")
## Analysis of Deviance Table ## ## Model: poisson, link: log ## ## Response: TOT.N ## ## Terms added sequentially (first to last) ## ## ## Df Deviance Resid. Df Resid. Dev Pr(>Chi) ## NULL 51 1071.4 ## D.PARK 1 680.55 50 390.9 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Intercept) D.PARK ## 4.3164849800 -0.0001058506
Y - Quantidade de veados infectados dentre os capturados
X - quantidade de vegetação aberta
plot(prop.veados ~ veados$openland,
ylim=c(0,1),
xlab = "Porcentagem de vegetação aberta",
ylab="Proporção de veados infectados")
nao.infectado <- veados$amostrado- veados$infectado y <- cbind(veados$infectado, nao.infectado) mod.veado <- glm(y ~ openland, data=veados, family="binomial") summary(mod.veado)
## ## Call: ## glm(formula = y ~ openland, family = "binomial", data = veados) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -3.8052 -1.1463 0.6707 1.3416 3.5483 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 2.4820 0.1736 14.29 <2e-16 *** ## openland -6.0336 0.4432 -13.61 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 376.18 on 31 degrees of freedom ## Residual deviance: 121.30 on 30 degrees of freedom ## AIC: 204.65 ## ## Number of Fisher Scoring iterations: 5
anova(mod.veado, test = "Chisq")
## Analysis of Deviance Table ## ## Model: binomial, link: logit ## ## Response: y ## ## Terms added sequentially (first to last) ## ## ## Df Deviance Resid. Df Resid. Dev Pr(>Chi) ## NULL 31 376.18 ## openland 1 254.88 30 121.30 < 2.2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Intercept) openland ## 2.482035 -6.033583
Teste de Shapiro-Wilks
shapiro.test(arv3$alt)
## ## Shapiro-Wilk normality test ## ## data: arv3$alt ## W = 0.96262, p-value = 0.3609
shapiro.test(arv$alt[arv$sp == "sp1"])
## ## Shapiro-Wilk normality test ## ## data: arv$alt[arv$sp == "sp1"] ## W = 0.9906, p-value = 0.9591
shapiro.test(arv$alt[arv$sp == "sp2"])
## ## Shapiro-Wilk normality test ## ## data: arv$alt[arv$sp == "sp2"] ## W = 0.97453, p-value = 0.3502
Teste de Bartlett:
bartlett.test(arv3$alt~arv3$sp)
## ## Bartlett test of homogeneity of variances ## ## data: arv3$alt by arv3$sp ## Bartlett's K-squared = 1.1109, df = 2, p-value = 0.5738
Teste de Levene:
leveneTest(arv3$alt~arv3$sp)
## Levene's Test for Homogeneity of Variance (center = median) ## Df F value Pr(>F) ## group 2 0.512 0.605 ## 27
RESĂDUO = Y(observado) - Y(dados ajustados)
A premissa de normalidade na verdade recai sobre os resĂduos do modelo.
resid(anova.sps)
## 1 6 11 16 21 26 ## 0.8395548 -0.5666556 3.1201253 -1.1287012 5.8661704 -12.1406501 ## 31 36 41 46 51 56 ## 6.7500566 -4.5005793 0.1541789 1.6065001 -2.6382600 -3.1208224 ## 61 66 71 76 81 86 ## 3.3224765 -3.6806771 -0.3024985 6.3022481 -1.3052922 7.4459310 ## 91 96 101 106 111 116 ## -2.4542276 -3.5688778 9.2780228 1.1454976 5.9847796 -9.3347719 ## 121 126 131 136 141 146 ## 4.4754638 -2.2630830 -5.5808309 -6.9952785 0.7849457 2.5052549
#residuals(anova.sps) # o mesmo que a função de cima
plot do modelo# regressĂŁo: alt ~ nutri plot(alt.nutr)
# ANCOVA modelo aditivo: alt ~ nutri + sps plot(alt.sp.nutr)
# ANCOVA modelo interativo: alt ~ nutri * sps plot(alt.sp.nutr2)
# GLM poisson: atrop ~ dist.park plot(mod)
# GLM binoimial: prop.infectado ~ openland plot(mod.veado)